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Abstract 

We introduce a class of exponential Runge-Kutta integration methods for kinetic equations. 
The methods are based on a decomposition of the collision operator into an equilibrium and 
a non equilibrium part and are exact for relaxation operators of BGK type. For Boltzmann 
type kinetic equations they work uniformly for a wide range of relaxation times and avoid the 
solution of nonlinear systems of equations even in stiff regimes. We give sufficient conditions in 
order that such methods are unconditionally asymptotically stable and asymptotic preserving. 
Such stability properties are essential to guarantee the correct asymptotic behavior for small 
relaxation times. The methods also offer favorable properties such as nonnegativity of the 
solution and entropy inequality. For this reason, as we will show, the methods are suitable 
both for deterministic as well as probabilistic numerical techniques. 

Keywords: Exponential integrators, Runge-Kutta methods, stiff equations, Boltzmann equa- 
tion, fluid limits, asymptotic preserving schemes. 

1 Introduction 

The numerical solution of the Boltzmann equation close to fluid regimes represents a major com- 
putational challenge in rarefied gas dynamics. In such regimes, in fact, the mean free path becomes 
very small and standard computational approaches lose their efficiency due to the necessity of us- 
ing very small time steps in deterministic schemes or, equivalently, a large number of collisions in 
probabilistic approaches. Several authors have tackled the problem in the past, and there is a large 
literature on the subject (see [2] |4j [TO] HB [H [36] and the references therein). Most standard tech- 
niques are based on domain-decomposition strategies and/or model reduction asymptotic methods. 
The direct time discretization of the Boltzmann equation, in fact, represents a challenge in such 
stiff regimes due to the high dimensionality and the nonlinearity of the collision operator which 
makes unpractical the use of implicit solvers. Exponential methods, like the so-called time relaxed 
discretizations |16j , combined with splitting strategies represent one possible way to overcome the 
problem. However, as discussed in various papers, the choice of the Maxwellian equilibrium trun- 
cation in the schemes was based more on physical considerations than on a direct mathematical 
derivation and it is an open problem to determine an optimal truncation criteria [SJ [32] ■ Despite 
this, these discretizations have been applied with success both in the context of spectral methods 
as well as Monte Carlo methods [31] . 

In this paper we propose a class of exponential integrators [22 [22] for the homogeneous Boltz- 
mann equation and related kinetic equations which are based on explicit Runge-Kutta methods. 
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The main advantage of the approach here proposed is that it works uniformly for a wide range of 
Knudsen numbers and avoids the solution of nonlinear systems of equations even in stiff regimes. 
Similarly to [16] , we derive sufficient conditions such that the resulting schemes can be represented 
as well as convex combinations of density functions including a Maxwellian term. This property is 
essential to achieve nonnegativity, physical conservations and entropy inequality. 

The starting point is the use of classical time spitting together with a decomposition of the 
gain term of the collision operator into an equilibrium and a non equilibrium part. Similar de- 
compositions for the distribution function has been used in [2] to derive unconditionally stable 
schemes and in [3D1 US EH H2] to develop hybrid Monte Carlo methods for kinetic equations. This 
decomposition of the Boltzmann integral has been introduced recently by Jin and Filbet [14] as a 
penalization method by a BGK-type relaxation operator to derive Implicit-Explicit Runge-Kutta 
schemes that overcome the stiffness of the full nonlinear collision operator. We recall that analo- 
gous penalization techniques based on linearized operators were previously used in the context of 
exponential methods for parabolic equations with applications to the Schrodinger equation [3TJ [25] 
and for kinetic equations [281 116) . In particular, in the present paper, we will show that for space 
homogeneous equations the Maxwellian truncation criteria introduced in [16] is equivalent to the 
penalization method in |14) . Let us mention that a study of the accuracy property of splitting 
methods in stiff regimes is beyond the scopes of the present paper. 

Even if we develop our schemes using the Boltzmann equation as a prototype model for its 
intrinsic difficulties and its relevance in applications, the methods applies to any large system of 
stiff ordinary differential equations of the form 



where e > is a small parameter, Y € M. N and the non-linear operator R(Y) is a dissipative 
relaxation operator as in [7] . Such operator is endowed with a nx N matrix Q of rank n < N such 
that QR(Y) = 0, V Y. This gives a vector of n conserved quantities y = QY. Solutions which 
belong to the kernel of the operator R(Y) — are uniquely determined by the conserved quantities 
Y = E(y) and characterize the manifold of local equilibria. Important examples of such dissipative 
relaxation operators arise in discrete kinetic theory, shallow water equations, granular gases, traffic 
flows and in general in finite difference/ volume discretizations of several kinetic equations. 
The methods here proposed use the following decomposition [T^] 



where N(Y) represents the non-dissipative non-linear part, A > is an estimate of the Jacobian 
of R evaluated at equilibrium and E(y) — Y is a simple dissipative linear relaxation operator. Note 
that, at variance with standard linearization techniques which operate on the short time scale, the 
operator is linearized on the asymptotically large time scale. 

This decomposition permits to apply exponential techniques which solve exactly the linear part 
and are explicit in the non-linear part. The use of such techniques, as we will see, is essential in order 
to achieve some unconditional stability properties of the numerical schemes [29] . Such properties 
are usually characterized as unconditional asymptotic stability and asymptotic preservation. 

Here we derive sufficient conditions for asymptotic stability and asymptotic preservation. This 
permits to introduce asymptotically stable and asymptotic preserving methods up to order 4 which 
are exact for BGK-type kinetic equations. 

The rest of the paper is organized as follows. First in Section 2 we present the Boltzmann 
equation and its fluid-dynamic limit. Operator splitting methods and various notions of stability 
are illustrated in Section 3. Next, in Section 4, we introduce the explicit exponential Runge-Kutta 
methods and derive conditions for asymptotic stability and asymptotic preservation. Examples of 



Y' = -R(Y), Y(t Q ) = Y Q 



(1) 



R(Y)=N(Y) + A(E(y)-Y), 



(2) 
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methods up to order 4 are also constructed. In Section 5 we describe numerical experiments for 
homogeneous problems and one application to the full Boltzmann equation. In a final section we 
discuss conclusions and perspectives of the methods proposed in this article. 



2 The Boltzmann equation 

We consider the Boltzmann equation of rarefied gas dynamics [B] 

d t f + v-V x f = Q{f,f) (3) 

with initial data 

f\t=o = Jo- (4) 

Here f(x, v, t) is a non negative function describing the time evolution of the distribution of particles 
with velocity v € M 3 and position x £ fl C M. d:c at time t > 0. In the sequel for notation simplicity 
we will omit the dependence of / from the independent variables x, v, t unless strictly necessary. 
The operator Q(f, /) which describes particles interactions has the form 

Q(/,/)= / BQv - v*\,n)[f(v')f(vi) -/(«)/(«.)] dv*dn (5) 

Jl 3 xS 2 

where 

u' = v + -(t> - v*) + -\v - v*\n, v'„ = v + -(v - v») - - \v - v*\n, (6) 

and B(\v — «*|, n) is a nonnegative collision kernel characterizing the microscopic details of the 
collision given by 

B(\v-v.\,n) = a( ^~ V *} ■ n) \v-V*\\ 
\\v-v*\ J 

with 76 [0,3). The case 7 = 1 is referred to as hard spheres case, whereas the simplified situation 
7 = 0, is referred to as Maxwell case. Note that in most applications the angle dependence is 
ignored and a is assumed constant [3]. 

The operator Q(f, f) is such that the the local conservation properties are satisfied 

mQ(f,f)dv=:(mQ(fJ))=0 (7) 
where m(v) = (l,i>, 4~) are ^ ne c °l ns i° n invariants. In addition it satisfies the entropy inequality 



2 

dt 



[ f log fdv= f Q(f,f)\ogfdv<0. (8) 

Jm.3 Jr3 



Integrating ([3]) against its invariants in the velocity space leads to the following set of non closed 
conservations laws 

d t (mf)+V x (vmf) =0. (9) 

Equilibrium functions for the operator Q(f, f) (i.e. solutions of Q(f, f) = 0) are local Maxwellian 
of the form 

M/(^,T) = ^-| FI exp(^f), (10 ) 

where p, u, T are the density, mean velocity and temperature of the gas in the x-position and at 
time t defined as 

(p,pu,E) = (mf), T = -^-(E — p\u\ 2 ). (11) 

3p 
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We will denote by 

U = (p, u, E), M[U] = M f , (12) 

clearly we have 

U=(mM[U}). (13) 

Now, when the mean free path between particles is very small the operator Q(f, f) is large and we 
can rescale the space and time variables in ([3]) as 

x' = ex, t' = et (14) 

to obtain 

d t f + v-V x f = -Q(f,f) (15) 

£ 

where e is a small parameter proportional to the mean free path and the primes have been omitted 
to keep notation simple. 

Passing to the limit for e — > leads to / — » M[U] and thus we obtain the closed hyperbolic 
system of compressible Euler equations for the macroscopic variables U 

d t U + V x F(U)=0 (16) 

with 

F(U) = (vmM[U}) = (pu, QU® u + pi, Eu + pu), p= pT, 

where / is the identity matrix. 

For small but non zero values of the Knudsen number, the evolution equation for the moments 
can be derived by the so-called Chapman-Enskog expansion [6] . This approach gives the compress- 
ible Navier-Stokes equations as a second order approximation with respect to e to the solution to 
the Boltzmann equation. 



2.1 Operator splitting and stability definitions 

Here we restrict ourselves to operator splitting based schemes. It is well-known, in fact, that 
most numerical methods for the Boltzmann equation are based on a splitting in time between free 
particle transport and collisions [3l [15] . Possible extension of the present theory to non-splitting 
schemes is actually under study and it will be considered elsewhere. 

Even if it is difficult to give a rigorous definition of asymptotic preserving scheme since the 
concept has been used for a long time in the physics and mathematics literature and may refer to 
different discretization parameters, here following |25 [ I26 [ 133 ] we formalize this notion for the time 
discretization of equation (j3|). 

Definition 1 A consistent time discretization method for of stepsize At is asymptotic pre- 
serving (AP) if independently of the initial data and of the stepsize At, in the limit e — > 
becomes a consistent time discretization method for the reduced system HI 6) ). 

Note that this definition does not imply that the scheme preserves the order of accuracy in t 
in the stiff limit e — > 0. In the latter case the scheme is said asymptotically accurate. 

As discussed above, the starting point in the solution of the kinetic equations is given by an 
operator splitting of ((3]) in a time interval [0, At] between relaxation 

dtf = -Q(f,f), (17) 

£ 

and free transport 

d t f + vV x f = 0. (18) 
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This situation is typical of Monte Carlo methods and of several other numerical codes used in 
realistic simulations. Even if this splitting, usually referred to as Lie- Trotter splitting, is limited 
to first order it permits to treat separately the hyperbolic free transport from the stiff relaxation 
step which often is of paramount importance in applications. 

Higher order splitting formulas can be derived in different ways (see |19j). Let us denote by 
7At(/) and C A t(f) the above transport and collision steps in a time interval [0, At] starting from 
the initial data fo, then the well-known second order Strang splitting |34j can be written as 

C A t/2(T At (C At/2 (f ))). (19) 

Unfortunately for splitting methods of order higher then two it can be shown that it is impossible 
to avoid negative time steps both in the transport as well as in the collision [19]. Higher order 
formulas which avoid negative time stepping can be obtained as suitable combination of splitting 
steps [11]. Note, however, that the appearance of negative coefficients or negative time steps in high 
order formulas may lead to several drawbacks in practical applications like the lack of positivity 
of the solution which makes very difficult their use in Monte Carlo schemes. 

As mentioned before the study of the accuracy properties of the different splitting methods in 
stiff regimes, although important, is beyond the scopes of the present manuscript and we refer to 
[241 135J for an error analysis of some splitting methods in stiff conditions. 

Now we can reformulate the asymptotic preserving property and prove that 

Proposition 1 A sufficient condition for a consistent time discretization method of stepsize At 
applied to the operator splitting approximation of given by {llfy - JTfy) . to be AP is that the time 
discretization of step Ji7] ). independently of the initial data and of the stepsize At, in the limit 
e — > projects the solution f over the local Maxwellian equilibrium M[Uq], Uq — (mfo). 

The proof of the above proposition is an immediate consequence of the fact that as e — > step 
(IT71) degenerates into the projection Ca*(/o) = M[Uq] which coupled with the transport step (TT51) 
originates a so-called kinetic approximation [8] to the Euler equation dl~6]) given by T A t(M [Uo]). 
We omit further details. 

In other words, Proposition Q] states that if the relaxation step ([17]) is AP then the whole 
splitting ([T7l) - (jl8[) is AP. Analogous results hold true for higher order splitting methods. 

In the sequel we will focus on the solution to the space homogeneous Boltzmann equation (fTTj) . 
In fact, most computational challenges related to the behavior of the full equation for small values 
of e depend on the time discretization of the homogeneous step. 

Of course AP is an important property in term of stability of the numerical scheme in stiff 
regimes. For implicit Runge-Kutta methods applied to (|17[) it has been shown in [33] that AP is 
equivalent to the notion of L-stability [H] . 

For general unconditionally stable schemes a weaker requirement is the notion of asymptotic 
stability (AS) [29]. Let us denote by f(t) and g(t) two solutions of ([P7]) corresponding respectively 
to the initial data fo and go such that Uq = (mfo) = (mgo). It can be proved that, for a suitable 
distance d(-, •), system ([17]) is contractive d(f(t),g(t)) < d(fo,go), and asymptotically stable since 
fit) — > M[Uo] and g{t) — > M[Uo] as t — > oo (see [6] for example) and so f(t) — g(t) —> as t —> oo. 
We refer to [37] for more details and examples of contractive metrics for problem ([17} in the case 
of Maxwell molecules. A particular metric is presented in Section [3~3l 

Let us denote by /" and g n , n > 1 the numerical solution at t — nAt obtained with a given 
time discretization method applied to (]17|) with initial data fo and go respectively. Now we can 
introduce the following definition. 

Definition 2 A time discretization method for ji?| ) is called unconditionally contractive with re- 
spect to the distance d(-, ■) if d(f , g 1 ) < e?(/o,.9o) holds for all fo, go such that (mfo) — (mgo) and 
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for all stepsizes At. Furthermore, it is called unconditionally asymptotically stable if f n — g n — > 
as n — ¥ oo independently of the step size At. 

Note that unlike contractivity, asymptotic stability is not related to a specific metric. Contrac- 
tivity of Runge-Kutta methods has been studied in [2Tj and it is well-known that such methods 
have limited order of accuracy. Clearly for implicit Runge-Kutta methods, asymptotic stability is 
closely related to ^-stability. 

For the sake of completeness we finally introduce the notion of entropic stability, namely schemes 
that preserve the entropy inequality jSJ. 

Definition 3 A time discretization method for |_?7| j is called unconditionally entropic if H{f n+l ) < 
H{f n ), where H(f) — J R3 flogfdv, independently of the step size At. 

The above monotonicity property for Runge-Kutta schemes is essentially equivalent to the so- 
called strong stability preserving (SSP) property which is often used when dealing with the time 
discretization of partial differential equations [T7J. Let us recall that Implicit Euler is the sole 
unconditional SSP method. All SSP Runge-Kutta methods of order greater than one suffer from 
some time-step restriction |20j . Thus, except for first order implicit Euler, the entropy inequality 
is not satisfied by high order unconditionally stable Runge-Kutta schemes applied to (fl7|) unless a 
suitable time step restriction is considered. 



3 Exponential methods 

Since we aim at developing unconditionally stable schemes, the most natural choice would be to 
use implicit solvers applied to (fTT)) . Unfortunately the use of fully implicit schemes for (fTT)) is 
unpracticable due to the prohibitive computational cost required by the solution of the very large 
non-linear algebraic system originated by the five fold integral appearing in Q{f,f) which has to 
be computed in each spatial cell at each time step in the inhomogeneous cases. We will see in this 
section a possible way to overcome these difficulties. 

First we rewrite the homogeneous equation (fTT)) in the form 

d t f=-(P(fJ)-nf), (20) 

£ 

where P(f, f) = Q(f, f) + [if and \i > is a constant such that P(f, f) > 0. Typically \i is an 
estimate of the largest spectrum of the loss part of the collision operator. Let us emphasize that 
most of the subsequent theory can be generalized to the case where P(f, f) is not strictly positive 
and fj, is an arbitrary nonnegative constant. However such an assumption makes the presentation 
easier and we will discuss possible generalizations later. 
By construction we have the following 

-(mP(f,f)) = (mf) = U. (21) 
M 

Thus P(f, /)//■< is a density function and we can consider the following decomposition 

P{f,f)/H = M[U] + g. (22) 

The function g represents the non equilibrium part of the distribution function and from the 
definition above it follows that g is in general non positive. Moreover since P(f,f)/fj, and M[U] 
have the same moments we have 

(mg) = 0. (23) 
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The homogeneous equation can be written in the form 

dtf = -9 + ~(M[U] - /) = H (?MlH - M[U]) + ^(M[U] - /). (24) 

The above system is equivalent to the penalization method for the collision operator recently 
introduced in [14] . Note that even if M[[7] is nonlinear in /, thanks to the conservation properties 
([7]), it remains constant during the relaxation process. The main feature of such formulation is that 
on the right hand side we have a stiff dissipative linear part /j,(M[U] — f)/e which characterizes the 
asymptotic behavior of / and a stiff non dissipative non linear part (P(/, /)//i — M[U])/e which 
describes the deviations of P(f, /)/ ' \i from M[U), or equivalently the deviations of the Boltzmann 
operator from a BGK-like relaxation term. 

We remark that the decision whether problem (flT|) . or equivalently (l20l) and ((24"]) . should be 
regarded as stiff or nonstiff does not depend only on the ratio fi/e but depends also on the chosen 
initial conditions. If the initial data is close to local equilibrium / = M[U] + 0(e), then the problem 
is clearly nonstiff. In fact, if / = M[f7] + efi we have 

Q(/, /) - e[Q{h, f) + Q(f, fx) + sQ{h,h% 

and so Q(f, f) = 0{e) and there is no need of using a specific stiff solver. For this reason, and the 
fact that the transport step (|18p drives the solution far from local equilibrium, in the sequel we 
concentrate our analysis to non equilibrium initial data. 

In such case the problem is stiff as a whole and a fully implicit method should be used in 
the numerical integration to avoid stability constraints of the type At = 0(e). On the other 
hand the linear part itself suffices to characterize the correct large time behavior of /. Therefore, 
instead of fully implicit methods, one should hopefully use methods which are implicit in the 
linear part and explicit in the non-linear part. One class of such methods is given by the IMEX 
Runge-Kutta schemes [TJ[33]. Note, however, that standard IMEX schemes may lose their stability 
properties since here also the explicit part is stiff. An alternative approach is based on the so-called 
exponential integrators where the exact solution of the linear part is used in the construction of 
the numerical methods EH [22] • 



3.1 Explicit exponential Runge-Kutta schemes 

In order to derive the methods it is useful to rewrite (|2~4"|) as 

W-M)e^ = l 
at e 

The above form is readily obtained if one multiplies (f2~4"f by the integrating factor exp(^i/e) 
and takes into account the fact that M does not depend of time. A class of explicit exponential 
Runge-Kutta schemes is then obtained by direct application of an explicit Runge-Kutta method 
to ([25]) . More in general we can consider the family of methods characterized by 

FW = e^/e F + Ag? g Aij{(iAt/£) (P^\ F<») _ M , 

(26) 



e — ' \ u 

3=1 



1 _ £ -c^At/e\ M n . = 1 



en+1 



= e-^r + ^±W i{ ,At/e) fg^g) -mA 

6 ti V i* J 

+ (l- e -^ At / £ ) M", 



(27) 
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where At is the time step, /" = f(t n ), M n = M(t n ), c t > 0, and the coefficients A i;j and the 
weights Wi are such that 

A ij (0) = a ij , Wi(0) = w i} i,j = l,...,v 

with coefficients a^- and weights Wi given by a standard explicit Runge-Kutta method called the 
underlying method. Various schemes come from the different choices of the underlying method. 
The most popular approaches are the integrating factor (IF) and the exponential time differencing 
(ETD) methods [2TJ [22J [29] . Since M n does not depend on time during the collision process in the 
sequel we will omit the index n. 

For the so-called Integrating Factor methods, which correspond to a direct application of the 
underlying method to (|25p . we have 

Aij(X) = aye~ (c *- Cj) \ i,j = l,...,v, i>j 

(28) 

Wi{X) = Wie-^-^, i = l,...,v, 

with A = fxAt/e. 

The first order IF scheme reads 

™ + l_„_£A, fn , fiAt^MM ( P(/ n ,/ n ) n f \ , f, ^tM 



which is based on explicit Euler. For such methods the order of accuracy is the same as the order 
of the underlying method. 

The Exponential Time Differencing methods are strictly connected with the integral represen- 
tation of (|2"5)) . In the general case the coefficients for ETD methods have the form 

Aij(X) = f e (1 - s ^ x Pl3 ( S )ds, i,j = l,...,v, i>j 
Jo 

Wi(\) = f eV-'^Pitfds, i = l,...,v, 
Jo 

where pi and pij are suitable polynomials. 

The standard first order ETD method based on explicit Euler in our case gives [21] 



/»+^=e- — / n + i^- tp\Z— j-X-LLJ., (30) 



where tp(z) = (1 — e z )/z. 



3.2 Time Relaxed methods 

A class of exponential methods for kinetic equations, the so-called time relaxed (TR) methods, has 
been introduced in [16] as a combination of an exponential expansion (or Wild sum) together with 
a suitable Maxwellian truncation. In this paragraph we show that these schemes included already 
decomposition (|24j) and can be derived directly from a suitable Taylor expansion of (|25|) . 
To show this, let us first introduce the change of variables 

r = 1 — exp(— jJit/e), 
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which, using the bilinearity of P(f, f), gives the equation 

1 



d_ 

a7 



(/ - M) 



1 -T 



= (P(/,/)-/iM) 



1 



M (l - r) 



2 ' 



(31) 



The application of an explicit Runge-Kutta scheme to (f3"Tj) with time step At = 1 — cxp(— /iAi/e) 
leads to a class of ETD methods. For example the first order scheme based on explicit Euler in 
the original variables yields 



/ 



ra+l 



7 n + 



£*Ai /^At\ fP{f n ,f n ) 



-<Pi 



M +(1- e" * )M, 



(32) 



where ^fc(z) = e~ z (l — e~ z ) k /z, k = 1, 2, 

Note that such scheme coincides with the first order exponential time relaxed method derived 
in [16] and differs from the standard ETD method based on explicit Euler. Higher order ETD 
methods can be derived as well simply applying higher order explicit Runge-Kutta methods to 
(l3~Tj) . Although interesting, here we do not explore further this class of schemes. 

Now let us consider a different approach by taking the Taylor expansion of (/ — M)/(l — r) 
around r = in (13T1). This leads to 



(/-M)/(l-r) = (/ -M)+r 



P(fo,fo) 



M 



P(P(foJo)Jo) + P(fo,P(fo,fo)) 



2M 



+ 0(r 3 ) 



where we have used the bilinearity of the operator P(f, /). 

If we compute all the terms in the expansion and use recursively the bilinearity of P(f, /) we 
can state the following 



Proposition 2 The solution to problem or equivalently V25\) or LSI]) can be represented as 

oo 

f(v,t) = (1 - r)f (v) + (1 - r)J2r k (.fj:(v) - M(v)) + tM(v), (33) 



k=l 

where the functions fk are given by the recurrence formula 

k 

/< 



f k+ i(v) = -i- -P(hJk-h)(v), k = 0, 1, . 

h=0 ^ 



(34) 



By truncating expansion (|33[) at the order m, and reverting to the old variables, we recover 
exactly the time relaxed schemes presented in |16) 

m 



fc=l 



which, using the fact that 



1 - e 



-fj,At/e 



Ed 



-fiAt/e\k 



)" = (!- e 



-fiAt/eyn+l 



k=0 
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can be rewritten in the usual form emphasizing their convexity properties 

m 

jn+l = e -»At/e - e -^ A */ £ ) fe /fc + (1 - e -^ At/£ ) m+1 M. (36) 

fc=0 

A remarkable feature of these methods is that the functions fk(v) are density functions with 
the same moments of the initial data (mfk) — (fo). Such property, together with unconditional 
nonnegativity and convexity of the weights, is enough to guarantee asymptotic preservation of the 
schemes as well as nonnegativity, contractivity and entropic stability (see [16] for details). 

Clearly TR schemes do not belong to the general family of methods defined by (j^5T) - (f?7)) . since 
they are based on the assumption that P(f, f) is a bilinear operator. 



3.3 Main properties 

In this section we will establish the main properties for a general exponential scheme in the form 
([26l) - (|27p . For some of the properties, like contractivity and entropic stability, we give proofs in 
the case of Maxwellian molecules. 

Since solutions to (|17[) are nonnegative and preserves the mass in our analysis we will restrict, 
without loss of generality, to probability density functions. 

Let us denote by T^R 3 ^ the class of all probability density functions / on R 3 , such that 



\v\ 2 dF(v) < oo. (37) 



We introduce a metric on ^(R 3 ) by 



1/(0-3(01 

j,y)= sup 

£GR 3 

where / is the Fourier transform of / 



d2{f,y) = sup (38) 



/(0 = e-**dF(v). (39) 



The metric (^(v) is nonexpanding with time along trajectories of the Boltzmann equation; 
that is, if / and g are two such solutions da(/(t), g(t)) < d 2 (fo,go)- The above property is a 
consequence of the fact that for Maxwell's molecules we have 

Q(f, /) = / <? f ' n) [/(«')/«) ~ /(«)/(«.)] dv * d n (40) 

where for any fixed unit vector e 

/ a(e ■ n) dn = S, 
Js 2 

with S > a constant. Taking fi = S and 

P(f, f)= [ o (^^T ■ n) f(v')f(vi)dv* dn, 

jRZxS 2 \\V-V*\ J 

we have Q(f, f) = P(J, f) + fif and [37] 

d 2 (P(f,f),P(g,g))<»d 2 (f,g). (41) 
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We refer to [37) for more information on the above metric and other contractive metrics for the 
Boltzmann equation in the case of Maxwellian molecules. 

Now let us denote by f n and g n the corresponding solutions obtained with an explicit expo- 
nential Runge-Kutta method. We have 



At 



F (i) _ Q (i) = e -c*^t/s(fn _ + _J2A lJ (LiAt/e)(P{F < -J\F^) - P(G U) ,G U) )) 



fn+i _ g n+i _ e -„At/etp _ 5 «) + At y> Wi(uAt/e)(P(F®,F®) - P(G®,G®)), 

e -f-f 

i—i 

and then 

da(f W,GW) < e^ At ^d 2 (f n ,g n ) + ^ £ \A t] {piAt/ e)\d 2 {F^\ G«) 

d 2 (/" +1 ,.9 ,l+1 ) < e-^ At / £ d 2 (r, 5 ") + ^Y]|m(/iAi/ e )|d2(i ?W ,G«). 

Let us denote by A = fiAt/e, A(X) the v x v matrix of elements |j4,j(A)|, w(X) the v x 1 vector of 
elements |Wi(A)|, J the v x 1 vector of elements d 2 {F^ l \G^) and e the vx 1 unit vector we can 
write 

(l-XA(X))d < E(\)d 2 (P,g n )e 
d 2 (f n+ \g n+1 ) < e- x d 2 (f n ,g n )+Xw(X) T d, 

where E(X) = diag(e" Cl \ . . . , e~ c " x ). 
Then we obtain 

d 2 (f n+ \g n+1 )<R(X)d 2 (P,g n ), (42) 

where 

R(X) = e- x + Xw(X) T (I-XA(X))~ 1 E(X)e (43) 
v-\ 

= e- A + ^A fc+1 u)(A) T A(A) fc ^(A)e. (44) 

fc=0 

To derive the last expression we expanded (/ — AA(A)) -1 by the geometric series and use the fact 
that A(X) is strictly lower triangular and so it is a nilpotent matrix of degree v. 
We can state [29 



Theorem 1 // an explicit exponential Runge-Kutta method in the form S26\) - f[27\ ) satisfies 

R(X) < 1, V A > 0, (45) 

with i?(A) given by \44ty then it is unconditionally contractive and unconditionally stable with respect 
to the metric d 2 (-,-). 

Note that for an IF method we have 

\Aij(X)\ < |a y |e-^-^\ |W,(A)| < K|e" (1 ~ Cl) \ 

thus we require 

= ci < c 2 . . . < c < 1, (46) 
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in order for the above quantities to be bounded independently of A. 
We have 



R{\) =e- x j l + ^V+V^e) . (47) 
Thus condition (|4"5)) is satisfied when 



fe=0 



v-l co , k+1 

1 + £A*+We< e * = ! + (48) 

k=0 k=0 V '' 

Now, if we consider an underlying Runge-Kutta method with a v x v non negative coefficient 
matrix A and a v x 1 non negative vector of weights u> then (|48|) holds if 

wTA " e "- (jfeTT)!' fc = - 1 '---^- 1 - 

The above condition is clearly satisfied if the underlying Runge-Kutta method is a j/-stage explicit 
Runge-Kutta method of order v. 
Thus we have proved 

Proposition 3 An explicit IF method is unconditionally contractive and asymptotically stable with 
respect to the metric •) if the underlying Runge-Kutta method is a v-stage explicit Runge-Kutta 
method of order v with nonnegative coefficients and weights satisfying \4&\) . 

As pointed out in [37J HH| examples of such methods are well-known up to v — 4 and the 
classical RK method of order four is the sole method with four stages. Moreover, it is also known 
that there does not exist explicit methods of order greater than four satisfying (|46[) with positive 
weights. For methods which are not of IF type there are no accuracy barrier, for example all time 
relaxed method satisfy immediately condition (|4"5)l . Other examples are reported in [2"2l |2"5] . 

We have 

Theorem 2 If an explicit exponential Runge-Kutta method in the form i26\) - [2l\ ) satisfies 

lim R(X) = 0, (49) 

A— >oo 

with i?(A) given by J^[ ) then it is asymptotic preserving. 

In fact taking go = M we have 

d 2 (f n+1 ,M) <R(X)d 2 (f n ,M), 

and so d2(f n+1 , M) goes to as A — > oo. 

For an IF method R(X) is given by (14"T[) and condition (14T)1) is always satisfied. Thus 

Proposition 4 An explicit IF method is asymptotic preserving if the underlying Runge-Kutta 
method satisfies J46| ). 

For practical applications it may be convenient to require that as A — > oo the numerical solution 
f n+1 and each level F^ of the IF method are projected towards the local Maxwellian without using 
explicitly the structure of P(f, /). It is straightforward to verify that this stronger AP property is 
satisfied if we replace condition (|4*51) by 

= ci < c 2 < . . . < c v < 1. (50) 
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We remark that the usual concept of stiff order [29] for exponential methods is in contrast with 
the latter strong AP property. For example for a stiff order one method the condition [32] 

1 - p- x 

implies that the Maxwellian term in (|27|i vanishes. So classical stiff order two ETD methods do 
not satisfy the AP property (see [29] for example). 

We conclude this section with a results concerning an important convexity property of the 
schemes. 

Theorem 3 // an explicit exponential Runge-Kutta method in the form i26\ )- [2l\ ) satisfies 

i — 1 — c ■ A 

J2 A ij( X ) ^ 'I ' > VA ^°> i=l,-.-,f (51) 

3=1 

X)Wi(A) < — VA>0, (52) 

luii/l Aij(A) > and Wi(X) > iften is unconditionally positive and entropic. 

In fact, the operator P(/, /) is nonnegative and convexity of the schemes is guaranteed if (I5T1) 
and Jig hold. Moreover since H(M) < H(f n ) and (see [38]) 

(Wl) < H(f), (53) 



by convexity we have also 

i-i 

< e- c<A i?(/ n ) + A^vly(A)F(FW) 

< #(/"), i = 1, . . . , i/ 

1/ 

H{f n+1 ) < e~ x H{f n ) + \J2Wi(X)H(FW) 

i=l 

+ (l-e- x -\J2w t (\))H(M) 



For convexity of an IF method we require 



2^aye CjA < , VA>0, i=l,...,u 



A 



u A i 

X^ eClA < < -^—i VA>0. 



A 

»=i 
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By Taylor expansion we obtain conditions 



3 = 1 

V 

i=l 

This allows to state the following 



X "<••'•/ < fe = 0,l,2,..., »=!,...,!/ (54) 



^iC- < fc = 0,1,2,..., (55) 



Proposition 5 ^4n explicit IF method is unconditionally positive and entropic if the underlying 
Runge-Kutta method has nonnegative coefficients and weights satisfying (CT )- f55P - 

Note that the above conditions on the choice of the underlying method are quite restrictive and 
that we are not using the bilinearity of P(f, f) which would lead to weaker constraints on aij and 
Wi. For example, if we consider the family of second order Runge-Kutta methods with two levels 
characterized by w± = 1 — w, ui2 — w and 021 = (2w) _1 with w € [0, 1], Proposition [5] is satisfied 
when 

w *-i>_+I fc = 0,l,... 
- 2 fe ' 

which implies WE [§, 1]. Examples of methods that satisfy convexity are the second order Midpoint 
or Runge method (w — 1) and the third order Heun method but not the classical fourth order 
Runge-Kutta method |18j . 

Remark 1 Convexity is an essential property if one wants to rely on Monte Carlo techniques for 
the computation of the collision operator. In fact, the resulting scheme is a convex combination of 
probability densities and then can be evaluated using the same methods described in [31]. 



3.4 Generalizations and implementation 

An essential aspect in the reformulation of the problem given by (j47|) is the choice of the value of 
the constant \i used in estimating the spectrum of the collision operator. Of course such constant 
can be chosen at each time step in order to improve our estimate. In the sequel we show different 
choices in the case of variable hard spheres. 

The choice of an upper bound for the loss part of the collision term leads to take fi = fi p where 

fip = sup / \v — v* | 7 /(v*) dv*. (56) 

Positivity is guaranteed since this choice implies clearly P(f, f) > 0. From a practical viewpoint 
computation of (1551) can be done at 0(N log N) for a deterministic method based on N parameters 
for representing f(v) on a regular mesh. This can be done using the FFT algorithm thanks to the 
convolution structure of the loss term in (|56[). Thus, in general, the computation of fi p will not 
affect the computational cost of the scheme. 

For Monte Carlo methods based on v\, . . . ,Vn particles one should estimate 

H P « max — ^ \v t - v 3 p . 

3 

To avoid the 0(N 2 ) cost it is a common choice to consider the following upper bound 

fL = 2 7 max|i>i - it| 7 > max -5- V* |^ - t>,| 7 , u=^-'S^v l . (57) 

i Vi i\ * — ' iv — * 

3 * 
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However, such positivity constraint on P(f, /) typically leads to overestimates of the true 
spectrum of the collision operator, especially in Monte Carlo simulations. A better estimate of \i 
would be given by the average collision frequency 



Ha = / \v - w*| 7 /(w)/(w*) dv« dv. (58) 

JR 3 JR 3 

This can be computed again at 0(N log N) cost in a deterministic setting whereas in a Monte 
Carlo simulation we have 

which, to avoid the quadratic cost, can be overestimated by 

Finally, as suggested in [2], /j, can be chosen as an estimate of the spectral radius of the 
linearized operator Q around the Maxwellian M. In fact 

Q{f, f) « Q(M, M) + Vg(M, M)(M - /) - VQ(M)(M - /), 

where VQ(M, M) is the Frechet derivative of Q evaluated at M. For example one can take 

Q(fJ) 



H s = sup 



f-M 



(60) 



Note however that the above estimate can be computed easily only in a deterministic framework. 

The choices fx = fi a or fi = ji s , although more accurate, pose the question of stability of the 
resulting scheme since they do not guarantee P(f,f) > 0. Note that, if we denote by P p (f,f ) = 
Q(fi /) + /%>/ an d byAp = [ipAt/e, for an arbitrary /i > 0, we can rewrite (|26|) - (|27|) as 



F« = e --Y"-(A P -A)2^(A)^)+A p g^(A)- Pp( ^ 0) ' FW) 



3=1 3 = 1 MP 



(61) 



l- e - c * A -A^^-(A) M, 



i =1 

? n+1 = e- A /"-(A P -A)^^(A)F^+A P ElF t (A) ^ lJ ' 

i=i i=i 

+ ( l-e~ A - A^Wi(A) ) M. 



(62) 



For stability now we must perform the same analysis of Section |3~31 on system (joTj) - (l6"2"j) . For brevity 
here we omit the resulting conditions which typically cannot be satisfied without introducing a 
stability restriction on the time step. 

To illustrate this let us consider the case v = 1 for IF methods. We have 

P ( f n f n \ 

f n+1 = (1 - A p + A)e- A /" + A p e~ A p{J ,J ' + (1 - e~ A - Ae- A ) M. (63) 

fJ-p 
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Convexity is guaranteed as soon as At satisfies 



At < , 

- I 1 

or equivalently if we take fj, > \i v — ej At which represents the lower bound for that makes the 
scheme unconditionally positive. On the other hand contractivity and asymptotic stability remain 
valid under the weaker restriction 

f + _L<i + e <-")>, P . 

Similar considerations hold for higher order IF schemes. In such cases the AP property is guaran- 
teed provided that the underlying method satisfy Proposition 2) 

Remark 2 Along this paragraph we have assumed fj, constant during the time stepping. In practice 
it is clear that the computation of [i v from i56\) on the initial data does not guarantee positivity of 
all terms in the vector f\ . In principle one can take fi = /i(t) and rewrite the exponential methods 
for a time dependent /j, from 

d( f - M)ei /o'cW* 1 i r * 

— £ = - (P /, /) - to K>) * (64) 

at e 

and then recompute /i p at each stage of the Runge-Kutta level. This procedure however may be 
quite expensive and in practice violation of positivity is rarely observed. In such circumstances one 
can set initially a numerical tolerance on /i p or simply repeat the computation with a larger value 
of n P . 



4 Numerical results 

In this section we perform some numerical tests for the exponential Runge Kutta schemes applied 
to the case of the full Boltzmann equation. We consider the first order IF scheme lp25)) and the 
second and third order IF schemes obtained using second order Midpoint and third order Heun 
methods respectively. All methods satisfy Proposition[5]and thus can be implemented using Monte 
Carlo strategies for the evaluation of the five fold collision integral [31] . Since we are interested in 
measuring the time discretization error of the different schemes we need to cancel the other sources 
of error in the computations. This is realized using a very large number of particles and statistical 
averages. 

We consider two different test cases: first the evolution of the fourth order moment in a space 
homogeneous case and then the heat flux behavior in a space inhomogeneous shock problem. 

4.1 Homogeneous relaxation 

As initial data we consider an equilibrium Maxwellian distribution with temperature T = 6, density 
q = 1 and mean velocity u = — 0.5. To this distribution we add a bump on the right tail along 
the x-axis. The bump is realized adding a Maxwellian with mass gi, = 0.5 g, mean velocity 
Uf, = 4 \jT and temperature T& = 0.5 T to the initial Maxwellian distribution. We consider the 
case of Maxwell interactions and hard spheres. The simulations are run till the equilibrium is 
approximately reached, which means t — 0.4 in the case of hard spheres and t — 0.8 for Maxwell 
molecules. The reference solution is computed by the Bird method [3] which converges toward the 
exact solution when the number of particles goes to infinity. 
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Figure 1: Li Error for the Fourth Moment Relaxation for the homogeneous relaxation problem 
with Maxwellian particles (left) and hard spheres (right). 
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In Figure [T] we show the Li error for the fourth order moment of the distribution function / 
for Maxwellian molecules and hard spheres. In each of the plots the error is depicted for different 
choices of the time step: respectively fiAt/e — 0.5, 1, 2 and 4. In the case of Maxwell molecules 
fj, = 1 while for hard spheres /j > 1 is an upper bound for the collision cross section computed 
using the simple choice (|57|) . As a consequence for the same values of fiAt/s the norm of 
the error is larger for hard spheres with respect to Maxwell molecules. Here we do not use any 
strategy to reduce the effect of possible overestimation of fj, as described in section 13.41 We leave 
this possibility to further research. The expected convergence rate is observed for all schemes. 

4.2 A shock wave computation 

We consider a Sod shock tube test with initial values 



The solution is computed using 150 grid points in space, the final time is t = 0.05. The transport 
step is solved exactly by particle transport as it is usual in Monte Carlo methods. As before we 
consider a very large number of particles and averaged the solution over several runs. Again the 
reference solution is obtained letting the time step go to zero and the number of particles to infinity 
using Bird's DSMC method. 

We report the heat flux profile for the first and the second order Runge Kutta integration factor 
scheme in figures [21 A first order splitting is employed for the first order IF while a second order 
Strang splitting is used in combination with the second order IF method. From top to bottom 
of the figures the Knudsen number values are e = 10 _3 ,5 10 -4 and 10 -4 . The time step is fixed 
equal to 10~ 3 . 

Both schemes show a good agreement with the reference solution for e = 10 -3 . Then, when 
the Knudsen number is halved we start to see some discrepancies between the profiles of the first 
order method and the reference solution. On the contrary the second order method is still in good 
agreement with the reference solution. When the Knudsen number becomes 1/10 of the time step 
both the schemes present larger errors but they are still able to catch correctly the departure from 
equilibrium. This is especially true for the second order scheme which is able to reproduce the 
correct profile with sufficient accuracy. 

5 Conclusions and further developments 

We have presented a general class of exponential schemes for the numerical solution of nonlinear 
kinetic equations in stiff regimes. The schemes generalize the class of schemes previously developed 
in [16j and share the fundamental property of asymptotic preservation. Even if the schemes have 
been developed in the case of the Boltzmann equation for dilute gases, extension of the schemes to 
other collisional kinetic equations which possess a smooth equilibrium state are straightforward. 
Let us also mention that decomposition (|22[) represents only one of the possible choices in order to 
linearize the collision operator close to equilibrium and then using it as a penalization factor in the 
construction of the numerical methods. For example a more accurate penalization can be obtained 
using the so called ES-BGK equilibrium function [23] , instead of the standard Maxwellian, which 
is well known to provide a better approximation of the collision dynamics. Let us finally mention 
that in principle the schemes can be extended to the non-splitting method case for ^ using the 
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Figure 2: Heat flux profile for first order (left) and second order (right) IF schemes. Top Kundsen 
number e = 10~ 3 , middle e — 5 10 -4 and bottom e = 10 -4 . At — 10~ 3 . 
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integral representation 

f(x,v,t) = /(x,v,0)e-^/ £ + ^ / e-^ t - s '>/ £ G(y(s))ds + ^ / e~> J - ( - t - s ^ e M(x,v,s)ds, 

£ Jo £ Jo 

where 

G(f,s) = --vV x f+^^--M. 
Here we do not explore further this direction and we leave it to future researches. 
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